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I. INTRODUCTION 

The nonlinear Schrodinger equation (NLS) models 
many kinds of wave phenomena. The NLS appears in 
diverse fields such as nonlinear optics , gravity waves 
on deep water 0, magneto-static spin waves 0, soli- 
tons in liquid crystals and magneto-sonic solitons in 
the atmospheric magneto- pause boundary layer It 
also describes the dynamics of the mean field of a weakly 
interacting atomic or molecular Bose-Einstein conden- 
sate (BEC) 0, 0, SJMj where it is known as the Gross- 
Pitaevskii equation 0, ^| . 

Here, we consider the steady state response of the mean 
field of a BEC to a potential step or a delta function po- 
tential, as modeled by the NLS. These potentials may be 
easily produced in present experiments on the BEC. The 
former may be realized by a detuned laser beam shined 
over a razor edge to make a sharp barrier, where the 
diffraction-limited fall-off of the laser intensity is smaller 
than the healing length of the condensate, so that the 
potential is effectively a step function. The latter models 
the response of the condensate to an impurity of a length 
scale smaller than the healing length, which could be re- 
alized by a tightly focused laser beam, by another spin 
state of the same atom, or by any other object, as for 
instance another alkali atom, confined in an optical trap. 
Moreover, the solution methods developed in this article 
may be generalized to arbitrarily complicated piccewise 
constant potentials. 

The complete set of stationary solutions to the NLS 
with a constant potential on the infinite line were dis- 
covered by Zakharov and Shabat 0, ^| . The station- 
ary solutions of the NLS under periodic and box bound- 
ary conditions has also been solved analytically |l4L Hoi , 
as well as the finite well Ha. The parabolic potential 
has been solved numerically [17|. The potential step has 
been examined theoretically and experimentally for the 
linear Schrodinger equation with a constant [lSL fl9| and 
oscillating |20j step. Klein examined similar problems 
experimentally by deflecting neutron beams with a vi- 
brating crystal [21|. In addition, symmetric steady-state 
solutions with a point-like impurity potential have been 
studied by Hakim and Taras-Semchuk |2^| in cer- 



tain limiting cases. Point-like impurity potentials have 
been studied extensively, such as with helium impuri- 
ties in a BEC l24[25l , BEC formation initiated by point 
like impurities |26( and impurity scattering in a BEC 
of sodium [2l}. The supcrfluid transmission of matter 
waves across various potentials has been studied [2^, . 
In addition, bound solutions to the NLS are the one- 
dimensional analog to the pinned vortex solutions which 
occur when a discontinuity is present in a two dimen- 
sional system, such as in a two dimensional high-T c su- 
perconducting system. 

In order to obtain the full set of stationary states in 
closed analytic form, we assume the BEC to be in the 
quasi-one-dimensional regime. When the transverse di- 
mensions of the BEC are on the order of its healing 
length, and its longitudinal dimension is much longer 
than its transverse ones, the ID limit of the 3D NLS 
is appropriate to describe the system [3(J. The ID NLS, 
with an external potential, V(x), may be written 
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d xx ^ + g\^\H + V{x)^, 



(1) 



where a harmonic oscillator confinement in the transverse 
directions with frequency lo has been assumed |3lj for 
atoms of mass M, the length has been rescaled accord- 
ing to units of the oscillator length, l/ lo = (fi/Muj) 1 / 2 , 
and energy rescaled according to units of the oscillator 
energy, Tiuj. The renormalized ID coupling, g = 2a s , 
where a s is the s-wave scattering length, characterizes the 
short-ranged pairwise interactions between atoms. The 
wave function or order parameter ^(x, t) has the physical 
meaning of ^{x, t) = y' p(x, t) exp[iip(x,t)], where p(x,t) 
is the longitudinal line density and the longitudinal su- 
pcrfluid velocity is given by v(x,t) = d(j)(x,t)/dx. Both 
attractive and repulsive atomic interactions, i.e., g > 
and g < 0, shall be considered. 

In the case where the harmonic oscillator length ap- 
proaches the s-wave scattering length, Iho — a s , the ID 
NLS no longer models the system and a one-dimensional 
field theory with the appropriate effective coupling con- 
stant must be considered instead |3ll |. Since a s is on the 
order of angstroms for typical BEC's, this regime is not 
relevant to the present study. Thus, it should be noted 
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that this study does not examine the Tonks-Girardeau 
regime |32j . where quantum fluctuations become impor- 
tant and the Gross-Pitaevskii equation no longer models 
the system. 

With the experimental demonstration of Feshbach res- 
onances in BEC's of dilute atomic gases [^3, 13- it is 
possible to alter the s-wave scattering length and, hence, 
the nonlincarity of the NLS. Near a Feshbach resonance, 
the scattering length becomes a function of a uniform 
background magnetic field. By altering the magnetic 
field, the effects of the nonlinearity can be experimen- 
tally controlled. It is therefore important to be able to 
characterize the complete set of solutions as a function 
of the interaction strength. 

This article is organized as follows. In Sec. |nj the full 
set of stationary solutions of the NLS along the infinite 
line are presented. In Sec. IIIII the stationary solutions 
influenced by a step potential are discussed. The effects 
of a point-like impurity are presented in Sec. IIVI Sec- 
tion shows the connection between the NLS solutions 
and the well-known solutions to the linear Schrodingcr 
equation. The results are discussed in Sec. IVII includ- 
ing a physical interpretation of all solution variables. In 
App. the special functions used throughout the arti- 
cle are briefly reviewed. Finally, in App. El the fact that 
the solutions discussed in Sec. [H] are indeed the full set 
of stationary solutions to the ID NLS with a constant 
potential [lj, [ljj is demonstrated formally. 



II. CONSTANT POTENTIAL 

A brief review is given of the stationary solutions of 
Eq. of the form, 



^(x, t) — R(x) exp[i(f)(x)] exp[— ipt] , 



(2) 



where p is the eigenvalue. There are several excellent 
references which consider stationary solutions in a con- 
stant potential [3(| HE OS |U |38| . Assuming a constant 
external potential V, the separation of Eq. (|TJ into real 
and imaginary parts gives 

Rd xx <t> + 2{d x R){d x <f>) = 0, (3) 
VR = fxR. (4) 



\d xx R+\{d x <f>) 2 R + gR 3 



(5) 



2 2 
Equation © can be integrated to give 

d x 4> = - , 
P 

where a is an undetermined constant of integration, and 
p(x) = R(x) 2 is the single-particle density. Substituting 
Eq. (jHJl into Eq. Q and integrating yields 

dp 
d.r 



2^/gp^-2(p-V)p 2 + Cp 



(6) 



Integrating a second time gives 
[ - 



: dp = 2x + x , (7) 



where C and xq are undetermined constants of integra- 
tion. It is shown in App. [B] that the only solutions to 
this integral equation are given by the Jacobian elliptic 
functions j3£||4(|. In App. |A"1 these special functions are 
briefly reviewed. The most general form of the solution 
is then given by 



p = Asn 2 (bx + x , k) + B . 



(8) 



where sn is one of the Jacobian elliptic functions, A is 
a density prefactor, b is a translational scaling, .To is a 
translational offest, k is the elliptic parameter, and B 
is a vertical density offset. The period of the density is 
given by 2K(k)/b, where K(k) is the complete elliptic 
integral of the first kind. Since the density is positive 
definite, the variables are restricted such that B > 
and B + A > 0. In Sec. IVII the relation between these 
variables and the mean number density, energy density, 
and momentum density are calculated and discussed for 
the nonlinear wave of Eq. © . 

It should be noted that since many of the solutions 
arc unbounded, the norm of the wavefunction remains 
unconstrained. This is in contrast to bound or localized 
solutions when the normalization, 



(p - p) dx = 1 



(9) 



may be used , where ~p is the average density. Alter- 
natively, it is possible to use a nonlinear scaling of the 
position and density, x — > ax and p — > a~ 2 p, to scale 
the maximum density on one side of the boundary to 
unity. However, throughout this paper the wavefunction 
remains unnormalized. 

One may then determine the variables p, a, and k, 
in terms of A, B, 6, and g by substituting Eq. (jSJ) into 
Eq. 10 and equating powers of the elliptic sn function to 
give 



1 



p 

, 2 



(b 2 + (A + W)g) + V . 



B(A + B){b 2 + Bg) 
A 



b 2 



9- 



(10) 

(11) 
(12) 



This leaves the eigenvalue, p, the constant of integration 
of the phase, a, and the elliptic parameter, k, determined 
up to A, B, b, and the interaction strength, g. Note 
that the fact that a enters into the equations only as 
a 2 implies that all nontrivial phase solutions, i.e., those 
for which a / 0, are doubly degenerate, as ±a lead to 
the same value of the eigenvalue, /i, without otherwise 
changing the form of the density or phase. We shall use 
the term trivial phase to refer to solutions for which the 
phase is spatially constant. 

In the following two sections these results are applied 
to pieccwisc constant potentials. In particular, the poten- 
tial step and the delta function potential arc examined. 
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III. POTENTIAL STEP 

In this section the complete set of solutions to the NLS 
with an external step potential of height Vb beginning at 

x = 0, 



V(x) = V o 0(x). 



(13) 



arc presented analytically, given the solution parameters 
on the negative x, or left, side of the step. In the fol- 
lowing two subsections, the general solution to the NLS 
with a potential step and some particular examples are 
discussed. 



A. General Solution 

Applying the two boundary conditions of continuity of 
the wave function and continuity of the derivative of the 
wave function gives the following five conditions, 



P (0+) 


= p(o~), 


(14) 


d xP (0+) 


= d xP (0-), 


(15) 


0(0+) 


= <j)(0~) + 2irn, 


(16) 


a(0+) 


= a(O-), 


(17) 


Mo+) 


= MO"), 


(18) 



where n is an integer. The first two conditions represent 
the continuity of the density, Eq. (|14fl . and the deriva- 
tive of the density, Eq. I|15|) . The next two conditions 
represent the continuity of the phase, Eq. (|16f) . and the 
derivative of the phase, Eq. ifTTjl . Equation (fT%|> demands 
that the eigenvalue be the same on either side of the 
boundary. Note that in Eq. (|10f> the eigenvalue on the 
right hand side is offset by Vb as compared to the left 
hand side. Since n only enters into the phase and does 
not effect such quantities as the eigenvalue and density, 
only the n = state is considered and, consequently, the 
phase, 4>, is continuous across the boundary. However, it 
is important to note that all solutions are of denumer- 
ably infinite degeneracy according to Eq. i|16|) . Since <fi 
is given by 



4>{x) 



a 



o P{x) 



dx + const. 



(19) 



continuity in the phase is easily achieved by setting the 
constant phase shift equal on either side of the boundary 
and, therefore, Eq. (jlGf) is satisfied. 

In the following derivation, it is assumed that the wave- 
function parameters on the left side of the step are known 
completely Therefore the density prefactor, Al, the ver- 
tical density offset, Bl, the translational scaling, b^, and 
the horizontal offset, xql are ah known, where the L sub- 
scribt refers to variables on the left side; an R subscript 
will refer to variables on the right side. In addition, the 
experimental parameters of the interaction strength, g, 
and the potential step height, Vb, arc both known. From 



the variables on the left and Eq. (|SJ), the density at the 
boundary, pl = p(0~), and its derivative at the bound- 
ary, d x pL = d x p(0~), can be determined. The eigen- 
value, p., and the phase constant, a, can be determined 
from Eqs. (JTUJ and lfTT|) . 

From Eqs. (jSJ and (|14|l the square of the Jacobian el- 
liptic sn function can be solved for, 



sn 2 (x 0R , 



PL 



B, 



A, 



(20) 



and, from Eqs. (|10|) and l|18|) . the horizontal scaling, b R , 



b R = 2(p-V )-(A R + 3B R )g. 



(21) 



These variables are substituted into Eqs. (|17|) and 11511 . 
using Eq. (f]~]~|) . to give 

a 2 L = B R (A R + B R ) 

x [2(p - Vb) - (A R + 2B R )g] , (22) 
{d xPL f = -4(S ft - p L )(A R + B R - p L ) 

x [2{p - Vq) - (A R + 2B R + p L )g] . (23) 

Equations 12 I'D and (|23|l are quadratic in A R and cubic in 
B R and can be solved analytically for A R and B R to give 
six solutions. The remaining variables on the right side 
can then be found by substituting the values of A R and 
B R into Eq. I|21|) . to find b R , and by taking the inverse 
Jacobi sn function of Eq. (|20(l to give 



xqr 



PL 



Br 



Ar 




(24) 



The full solution is then completely known. It is there- 
fore possible to completely describe the system analyti- 
cally knowing only the parameters on one side of the step. 
This not only introduces computational ease in evaluat- 
ing solutions, it also provides all possible solutions, most 
of which cannot be determined using purely numerical 
methods. In the following section, specific examples of a 
BEC in the presence of a step potential will be examined. 



B. Particular Examples 

The solutions to the potential step problem can be di- 
vided into two categories. The eigenvalue, p, can be large 
enough that particles are free to move across the bound- 
ary and the nonlinear analog to a transmitted wave for 
the linear Schrddinger equation becomes possible. How- 
ever, if the eigenvalue is too small, then the wavefunction 
must decay under the step. 

When p is larger than the effective potential, 



v eff(x) =V + gp{x) , 



(25) 



the wave can be transmitted across the boundary. Note, 
in the case of an attractive interaction, g < 0, the eigen- 
value can be less than the step height, Vb, and for a re- 
pulsive interaction, g > 0, the eigenvalue must be strictly 
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FIG. 1: Stationary solutions to the NLS with a potential step 
of form V(x) = Vo9(x). These solutions, which are the non- 
linear analogs of the continuum of linear stationary scattering 
states, exhibit a large deviation from the traditional linear so- 
lutions. Shown are particular examples of (a) the density and 
(b) the phase for a repulsive interaction strength and (c) the 
density and (d) the phase for attractive interaction strength. 



FIG. 2: Stationary solutions to the NLS with a potential step 
for which the wavefunctions decay, which are the nonlinear 
analogs to evanescent waves. Shown are particular examples 
of the densities of nonlinear waves with (a) repulsive interac- 
tion strength, (b) attractive interaction strength, and (c) and 
(d) with an attractive interaction strength that decay on the 
lower side of the step. 



greater than the step height. In Fig.^a) the density of 
a nonlinear state with a repulsive interaction strength is 
shown. A step of height Vq = 1, positioned at x = 0, 
and a condensate with an interaction strength of g = 0.2 
and eigenvalue of // = 2.404 were used. Notice that 
the increased interaction strength and nonlincarity has 
caused the peaks of the wavefunction to become much 
broader than in the linear case. The phase that corre- 
sponds to this density is shown in Fig. ^b). Figure ^c) 
shows a similar solution but with an attractive interac- 
tion strength. This potential is again given by a step 
with height of Vq = 1, positioned at x = 0. An interac- 
tionstrength of g = —0.2 and eigenvalue of /i = 0.98 was 
used. In this case, the peaks have instead narrowed due 
to the attractive interaction. The phase that corresponds 
to this density is shown in Fig. ^d). 

When the eigenvalue is less than the effective potential, 
fj, < V e ff(x), the wavefunction must decay under the 
step. In Fig.[2Ia) the density of a nonlinear solution with 
a repulsive interaction strength that decays as it crosses 
the boundary of the step is shown. A step with height of 
Vq = 1, positioned at x = 0, and a condensate with an 
interaction strength of g = 0.12 and eigenvalue of /i = 0.5 
were used. Figure |2Ib) shows a similar solution but with 
an attractive interaction strength. This potential is again 
given by a step with height of Vq = 1, positioned at 
x = 0. An interaction strength of g = —10 and eigenvalue 
of yu = —49 were used. For both wave functions, the 
phase is necessarily trivial, since all wavefunctions that 
approach zero at infinity must have A = — B, and hence 
from Eqs. © and (|llf) . the phase is constant. In order 
to stress the importance of the effective potential, and 
not just the step potential, Figs. EJc) and (d) show the 



density of two nonlinear wavefunctions with attractive 
interactions, g = —1, that decay on the lower side of 
the potential. Both wavefunctions have an eigenvalue 
of /.t = —0.5 and a potential of height of Vb = 0.01 and 
Vq = 1 were used for Fig. (J2tc) and Fig.[5Id), respectively. 

Thus the NLS with a potential step has solutions to 
the step potential that provide a connection between the 
linear solutions and a wide variety of exotic nonlinear 
wavefunctions, as shall be discussed in Sec. |Vj 



IV. POINT-LIKE IMPURITY 

In this section, the case of a potential of form 

V(x) = V 5(x) (26) 

is considered. Such a potential models an impurity which 
deforms the constant background potential on a length 
scale much less than that of the healing length. Given 
the state of the system on the negative x, or left side 
of the impurity, the variables on the right side of the 
impurity are determined. A negative value of Vo repre- 
sents an attractive impurity ; such as due to defects in 
hydrogen-bonded chains (4ll l42fl , and a positive value of 
Vq represents a repulsive impurity, such as with helium 
atoms in a BEC [mI^. 



A. General Solution 

The boundary conditions for an impurity are sim- 
ilar to those for the potential step, except that the 
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derivative of the wavefunction experiences a disconti- 
nuity at the boundary. Therefore, it is necessary that 
Eqs. (HJj, CHI), and (32} must still be satisfied, as 
well as 

d xP (o + ) - a x p(o-) = -4 P (o)y • (27) 

It is again assumed that all variables on the left side 
of the impurity are known as well as the experimen- 
tal parameters of interaction strength, and impurity 
strength, Vq. Using a treatment that is exactly anal- 
ogous to that for the step function, all of the parame- 
ters on the right side, given those on the left, are an- 
alytically determined. The only difference is that in 
Eqs. (J2TJ, (22}, (2j|, and (22}, the quantity d x p L must 
be replaced with (d x pL — iVopL,) an d {p — V ) must be 
replaced with p. It is therefore possible to completely de- 
scribe the system analytically knowing only the parame- 
ters on one side of the impurity. In the following section, 
examples of a wavefunction subject to an impurity are 
examined. 



B. Particular Examples 

For the delta function potential, both symmetric and 
nonsymmctric wavefunctions are possible. Of particular 
interest are the symmetric wavefunctions in the k = 1 
limit of the Jacobian elliptic functions. In this case all 
solutions become hyperbolic trigonometric functions with 
a localized change in the density around the impurity and 
no oscillations at ±00. Solutions of this type we term 
localized 0|. Due to the form of the solutions, there are 
four different solution types. The possible wave functions 
are then 

p = a 2 + (1 - a 2 ) tank 2 {\ZT~c7 I \x\ + x ) , (28) 

p = a 2 + {1 - a 2 ) coth 2 (y/T~c^\x\ + x ) , (29) 
p = b 2 sech 2 (b\x\ +x ), (30) 
p = b 2 csch 2 (b\x\ + x ) , (31) 

where the transnational offset, xq, is determined by the 
impurity strength, Vq, and the density has been normal- 
ized according to Eq. @. Equations (|28J) and (|29|l arc 
valid for repulsive interactions, while Eqs. (|3l)|l and 131|) 
are valid for attractive interactions. 

Figure[3]shows the possible wavefunctions for repulsive 
interactions. An example of the solution described by 
Eq. (|28|l is plotted in Fig. Ola) with a repulsive impurity 
of strength Vo = 0.5. This may be interpreted as a single 
dark soliton bound by an imp urity. This is similar to the 
solution found by Hakim [22| for a soliton that is moving 
with an impurity. In addition, a bound state of two dark 
solitons, Fig. |3Jc), can be created when the strength of 
the impurity is attractive and exactly balances the re- 
pulsion between the two dark solitons. Figure|2Ic) shows 
the hyperbolic cotangent function solution with an im- 
purity strength of Vq = —0.5; this may be interpreted as 



a deformation of the ground state constant solution to 
the NLS with a constant potential. In all plots an inter- 
action strength of g = 1 and phase constant of a = 0.5 
were used. 

It should be noted that there exists a bound state of 
a repulsive condensate with an attractive impurity. This 
solution is given by, 

-csch^ta + cotlrV^)), (32) 
9 b 

where the interaction strength, g has specifically been 
included and b must be determined such that the density 
is normalized to unity. These requirements place a limit 
on how repulsive the interaction may become and is given 
by, 

ffraax = -4Vo , (33) 

where g max is the most repulsive interaction the conden- 
sate may have. If the interaction is increased past this 
point, the condensate will spill away from the impurity 
and will no longer be bound. 

The set of symmetric localized solutions for the case of 
attractive interactions do not allow for nontrivial phases, 
in contrast to the case of repulsive interactions. For at- 
tractive interactions, g < 0, the hyperbolic secant func- 
tion solution, Eq. (201, is valid for both V > 0, Fig.QJa), 
and Vq < 0, Fig. Bib), where potential strengths of 
Vo = 0.9 and Vq = —0.9 were used, respectively. These 
solutions may be interpreted single bright soliton, 
which is the ground state solution to the 1D-NLS, de- 
formed by an impurity. The hyperbolic cosecant function 
solution, Eq. I|31[) . is only valid if Vo < and is similar in 
form to the hyperbolic secant solution of Fig. EJb) • The 
two solutions types are degenerate for Vo < 0, with an 
eigenvalue of p = —b 2 /2. In Figs, ^a) and (b) an in- 
teraction strength of g = — 1 and translational scaling of 
b = 1 were used. 

Nonsymmctric wavefunctions which oscillate at infin- 
ity are also possible and come in two forms: oscillations 
on one side of the delta function and oscillations on both 
sides. Figure [5] shows two possible nonsymmetric wave- 
functions subject to a delta function positioned at x = 
that oscillate on one side of the potential. The density 
of a wavefunction with a repulsive interaction strength of 
g = 0.21 and eigenvalue of p = 2.4, distorted by a delta 
function, Vq = 2, is shown in Fig.[S{a), where the left side 
reproduces the hyperbolic tangent function of Eq. I|28|l . 
In Fig. [SJb), an attractively interacting, g = —50, wave- 
function with an eigenvalue of p = —50, distorted by a 
delta function, Vo = 10, is shown that appears similar 
to the evanescent wavefunctions of Fig. [2J that decay be- 
neath a step. Note that Fig. |S{a) has a nontrivial phase 
while the phase of Fig. EJb) is trivial. 

Figure [B] shows two possible nonsymmetric wavefunc- 
tions subject to a delta function positioned at x = with 
strength Vq = 2 that oscillate on both sides of the poten- 
tial. A repulsive interaction strength produces the char- 
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FIG. 3: Localized, symmetric solutions to the NLS with re- 
pulsive interaction strength in the presence of an impurity, 
V(x) — VqS(x). Shown are particular examples of (a) the 
density and (b) the phase of a dark soliton bound by a re- 
pulsive impurity, (c) the density and (d) the phase of a pair 
of dark solitons bound by an attractive impurity, and (e) the 
density and (f) the phase of a supercurrent deformed by an 
attractive impurity. Note that (a) and (b) may also be inter- 
preted as deformations of a supercurrent. 
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FIG. 4: Localized, symmetric solutions to the NLS with at- 
tractive interaction strength in the presence of an impurity. 
Shown are particular examples of (a) the density of a bright 
soliton, which is the ground state solution to the NLS, de- 
formed by a repulsive impurity and (b) the density of a bright 
soliton deformed by an attractive impurity. 



acteristic widening of the pulse peaks, Fig.UJa). The cor- 
responding phase is given in Fig. [fJJb). An attractive in- 
teraction strength creates a narrowing of the pulse peaks, 
Fig-EJc)- The corresponding phase is given in Fig.EJd). 
The condensates in Figs. EJa) and (c) are characterized 
by eigenvalues of /i = 2.4 and fx = —1.3, respectively. 




J -5 5 

Position (units of l ho ) 



5 0.5 

Position (units of l h0 ) 



FIG. 5: Solutions to the NLS with an impurity that oscillate 
on one side of the potential. Shown are particular examples 
of (a) the density of a nonlinear wave with repulsive inter- 
action strength and (b) the density of a nonlinear wave with 
attractive interaction strength. These waves have no analog 
with the solutions of the linear Schrodinger equation. 
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FIG. 6: Nonsymmetric solutions to the NLS with an impurity. 
Shown are particular examples of (a) the density and (b) the 
phase of a nonlinear wave with repulsive interaction strength 
and (c) the density and (d) the phase of a nonlinear wave with 
attractive interaction strength. These solutions are the non- 
linear analogs to the continuum of linear stationary scattering 
states. 



V. LINEAR LIMITS 

In this section it is shown how the solutions to the 
NLS connect to the solutions of the linear Schrodinger 
equation. There are two distinct types of possible waves. 
If the wave has enough energy it is possible to make the 
wave propagate through space. However, if it does not 
have enough energy, the wave can carry no current and 
is refered to as an evanescent wave. These waves must 
decay. In the next two sections, these two types of waves 
are discussed concerning their role as the linear limit of 
the nonlinear solutions. 



A. Transmitted Waves 

In a linear system, if the energy of the system is greater 
than the potential energy, then the wave can be transmit- 



7 



ted through space as sine waves. The nonlinear Jacobian 
elliptic sn function waves are the nonlinear analogue of 
the sine waves of a linear system. In this section, it is 
shown how the nonlinear solutions, in the linear limit, 
recreate the linear solutions of propagating waves. 

For the linear case, the usual representation of the sys- 
tem by incident, transmitted, and reflected waves is given 
by, 



^ik L x _ 
rp ik R x 



— ifit 



(34) 
(35) 



where k L = y / 2(p -V L ), k R = y/2(n - V R ), pt is the 
eigenvalue of the Schrodinger equation, R is the reflection 
coefficient, T is the transmission coefficient, and the inci- 
dent wave is assumed to be coming in from the left. The 
potentials, Vl and Vr, and the scattering coefficients, R 
and T, are determined by the type of boundary. For 
the case of a potential step, the potentials are given by 
Vl = and Vr — Vq. For a delta function potential, the 
potentials are given by Vl = Vr = 0. The wavefunction 
can be alternately described by an amplitude and phase 
as follows: 



d x (j) 



PL = (l + r) 2 -4rsin 2 (V^Ma;-s/2), (36) 
(r-l)(r + l)V27Z 



(1 + r) 2 -4r sin 2 ( v / 27ta; - s/2) 



, (37) 



where 



*z,(x,t) = s/p L (x) exp^O^e-^*, (38) 
R = re ls , (39) 

with r and s real. The nonlinear solution is then con- 
nected to the linear solution by 



x a = -s/2, 

b 2 = 2/i, 

A = -4r, 

B = (r + 1) 2 . 

The transmitted wave can easily be given by 
PR = t 2 , 



d x <t>i 



2fi 



where 



$R(x,t) = \/pr(x) exp[i<f) R (x)]e v * 
T = te iw , 



(40) 
(41) 
(42) 
(43) 

(44) 
(45) 



(46) 
(47) 



with t and w real. The nonlinear solution is then con- 
nected to the linear solution by 



x = 0, 

b 2 = 2fi. 

A = 0, 

B = t 2 . 



(48) 
(49) 
(50) 
(51) 



The nonlinear solutions when the eigenvalue is greater 
than the effective potential are therefore adiabatically 
connected to the linear transmitted wave solutions. The 
next section connects the decaying evanescent waves with 
nonlinear solutions. 



B. 



Evanescent Waves 



If the eigenvalue, /x, is less than the effective potential 
then the wave function must decay In the linear case, 
the decay is precisely exponential. The Jacobian elliptic 
sn function solution, Eq. JHJ), of the NLS can provide the 
appropriate exponential decay. In the limits g — > and 
B, xq — > +oo under the constraints A=— B and A g/b 2 = 
1, Eq. ijHJ) gives 



lim B sech (bx 

B ,xq — >+oc 



xa) = Pq e 



- bx 



(52) 



where po is the density at the boundary and is given by, 



Po 



e 



4B 

2x 



(53) 



When the limit that g approach zero is not enforced, 



p(x) oc 



1 



(N+ exp+ bx +N- exp-^) 2 ' 



(54) 



where N + and N_ are constants related to the magni- 
tude and sign of the interaction strength. The decay of 
the density is, therefore, not strictly an exponential de- 
cay. It is interesting to note that all decaying solutions, 
whether linear or nonlinear, must possess a trivial phase. 
This is due to the restriction that A=—B, and, hence 
from Eqs. (|TT|> and ©, the phase must vanish. This 
is also consistent with the physical interpretation of the 
phase since a nontrivial phase corresponds to a supcrfluid 
velocity and the velocity must vanish if the wave cannot 
be transmitted. 



VI. DISCUSSION AND CONCLUSIONS 

The full set of stationary states of the mean field of 
a Bose-Einstein condensate, modeled by the nonlinear 
Schrodinger equation in one dimension, in the presence 
of a potential step or pointlike impurity were presented 
in closed analytic form. Non-decaying solutions were di- 
vided into two categories: localized soliton-like solutions, 
and solutions that oscillate out to infinity. The localized 
solutions are of a purely nonlinear character, as they have 
no linear analog. The oscillating solutions, on the other 
hand, were shown to be adiabatically connected to the 
solutions to the linear Schrodinger equation. 

The localized solutions present novel wavefunctions. 
With a delta function potential, the localized solution can 
be interpreted as a single bright or dark soliton trapped 



8 



by the impurity. In addition, it was shown that an im- 
purity can also bind a soliton pair. If the impurity is at- 
tractive, the natural repulsion between two dark solitons 
can be exactly canceled by the attraction of the impu- 
rity, while if the impurity is repulsive, it can balance the 
natural attraction of in-phase bright solitons. Since these 
solutions are conjected to be stable (see below), they are 
excellent candidates for the experimental realization of 
stationary excited states of a Bose-Einstein condensate. 
In addition, the maximum repulsive interaction strength 
of the condensate with an attractive impurity that allows 
for a bound state has been determined. 

The oscillating solutions to the NLS, despite being 
adiabatically connected to oscillating solutions to the 
Schrodinger equation, have very different properties due 
to the concept of an effective potential. For the attrac- 
tive interaction solution for an evanescent wave decaying 
under a step, as illustrated in Fig.^b), the eigenvalue is 
larger than the effective potential in the regions of high 
density and is less than the effective potential in regions 
of low density. Figures EJc) and (d) show more radical 
deviations from the linear solutions since the wavefunc- 
tions decay on the lower side of the step. 

It is possible to characterize the general solution to the 
NLS, Eq. |(HJ), in terms of physical quantities such as the 
mean linear number density, n , mean energy density, £ , 
and mean momentum density, V. The densities are given 
by 



= B + A[ —t - 



/J_ E{k) x 
U 2 k 2 K{k)J 



(55) 



£ 
V 



3B 2 g 
nfx H 



a . 



b 2 A 



\(li-V){n-B), (56) 
(57) 



where fi, a, and k are given by Eqs. I|1(J|) . i|llfl - an d tUt ; 

respectively. In general, the number, energy, and mo- 
mentum densities must be calculated separately for the 
left and right sides of a boundary. These densities can be 
used to determine the variables A, B, and b of Eq. ||SJ), 
leaving only the translational offset, xq, as a free variable 
determined by the boundary conditions. It should be 
noted that the factor multiplying A in the mean number 
density, Eq. I|55|l . approaches one half when k approaches 
zero and approaches one when k approaches unity. This 
is to be expected since the mean number density of a 
linear wave is given by B + A/2 and the mean number 
density of an extremely nonlinear wave is given by B + A. 
The mean energy density can easily be calculated from 



c — y 2 
£ = n\i- -p 2 , 



(58) 



and so the second and third terms on the right side of 
Eq. I|56|) are due to this nonlinear correction. Since the 
mean momentum density, Eq. I|57|) . is equal to a, the 
mean momentum density must be equal on both sides 
of the boundary due to the boundary condition on a, 
Eq. (|17|l . In addition, the momentum density as a func- 



tion of position is also given by a and so the momentum 
density is equal everywhere. 

The healing length, £, of the NLS in the quasi-one- 
dimensional regime, where the transverse dimensions are 
trapped by a harmonic potential of frequency u>, is given 
by 



e 



/2 



(59) 



where n is given by Eq. I|55|) . Since the mean number 
density can vary across the boundary, it is possible for 
the condensate to have a different healing length on ci- 
ther side of a boundary. Since the speed of sound in 
the condensate is inversely proportional to the healing 
length, the speed with which phonon-like excitations can 
travel vary as they cross the boundary. 

While finding the complete set of solutions to Eq. 
with an impurity or step potential provides much infor- 
mation about the system, only stable solutions are exper- 
imentally observable. Previous works have examined the 
stability of stationary states for a constant external po- 
tential (see, for example, [M El El El El El El El ) , 
as well as for periodic and harmonic potentials (see, for 
example, El Ell)- Most studies are ultimately numeri- 
cal: linear stability can be solved in a few special cases, 
while nonlinear stability is analytically intractable. For 
a constant external potential, single bright and dark soli- 
tons and dark soliton trains are stable. A finite num- 
ber of bright solitons may form bound states, as for ex- 
ample order n solitons (n > 1). Bright soliton trains 
are always unstable, but may be experimentally stable, 
in that their lifetime is much longer than experiments, 
which typically require stability timescales of from mil- 
liseconds to seconds. Bright soliton trains which have a 
phase difference A(f> between adjacent peaks such that 
— 7r/2 < A<f> < 3ir/2 exhibit this experimental stabil- 
ity, with the lifetime being longer the closer A<f> is to ir. 
Bright soliton trains with — n/2 < A<f> < ir/2 arc unsta- 
ble but become quasi-periodic in time in a finite system. 

Based on these known results from the case of a con- 
stant potential, as well as the from the stability analysis 
with an impurity performed by Bogdan et al. E3j the 
stability of an attractive condensate with an impurity is 
as follows. According to Bogdan, the bound state of two 
bright solitons, as in Fig. 0Jb), are stable, since, so long 
as they are strongly overlapping, they will be in phase 
(A(f> = 0) and remain bound to the impurity. This is 
also the ground state of the system. However, the kind 
of solution shown in Fig.^Ja) is unstable Ej|. We then 
conject on the stability of a repulsive condensate that 
is not bound and whose density approaches a nonzero 
constant at infinity. Localized solutions in the case of 
repulsive nonlinearity are obviously stable in the cases of 
Fig-EI a )-(b) and|!Je)-(f), since they are the ground state. 
The bound pair of dark solitons, illustrated in Fig.^c)- 
(d), should be likewise stable, so long as the impurity is 
sufficiently strong. All of these solution types, except for 
attractive solitons bound by a repulsive imputiry, are ex- 
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pected to be experimentally observable in finite systems, 
such as an elongated harmonic trap. An excellent anal- 
ysis of the stability of solitons pinned with impurities is 
given by Bogdan, et. al. . 

The stability of soliton trains is a less certain issue. In 
the repulsive case, the central question is whether or not 
the phase locking of the individual solitons in the train is 
destroyed by the impurity or potential step; if it is, they 
may become unstable in the region of the discontinuity 
in the potential. In the case of bright soliton trains, un- 
less there is a strong phase difference between the peaks, 
they will attract and become unstable; otherwise, the 
discontinuity should not present a source of instability, 
since bright solitons adjust themselves to perturbation 
by emission of a small fraction of the total wavefunc- 
tion [4^. In order to perform numerical studies, the so- 
lutions would have to be quantized on a ring, in order to 
provide a finite domain for simulation. Such a stability 
study presents a subject for future research. 

We emphasize that neither the idea of left and right 
traveling waves nor that of reflected plus incident waves 
apply to nonlinear wave equations. This is important 
since one cannot create wavepackets from linear combina- 
tions of these solutions. Instead, these solutions already 
contain the wavcpackct-like solutions, or solitons, that 
are necessary to describe the system; moreover, solitons, 
unlike wavepackets, are nondispersive. Time dependent 
nonlinear scattering remains an open question that can 
certainly be addressed via numerical studies. In general, 
the stationary solutions to the NLS give physical insight 
into its dynamics, without which numerical solutions may 
be difficult to interpret. Perhaps more importantly, us- 
ing the general nature of the solutions to the cases of a 
step function and an impurity, it is possible to describe 
all stationary states to piecewise constant potentials. 

In conclusion, we have analytically solved for all sta- 
tionary solutions to the nonlinear Schrodingcr equation 
with a delta function or a step function potential. This 
models the steady state behavior of the mean field of a 
Bose-Einstein condensate in the presence of an impurity, 
or of a potential step created by, for instance, a laser pass- 
ing over the edge of a razor blade. Novel wavefunctions 
were found, including solitons trapped by the impurity 
and the nonlinear analog of transmitted and evanescent 
waves. 
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TABLE I: Limits of the Jacobian elliptical functions and in- 
tegrals l4?3l . 





u n 

n — u 


U 1 

ft, — 1 


sn(i/, k) 


sin(ii) 


tanh(it) 


cn(ii, fc) 


cos(ii) 


sech(ti) 


dn(u, k) 


1 


sech(u) 


ns(w, k) 


csc(w) 


coth(u) 


ds(w, k) 


csc(m) 


csch(u) 


cs(u, k) 


cot(u) 


csch(u) 


K(k) 


tv/2 


00 


E(fc) 


vr/2 


1 



APPENDIX A: JACOBIAN ELLIPTIC 
FUNCTIONS 

A brief review of the Jacobian elliptic functions 0, |5(j 
is given. Of the 12 elliptic functions, there are only six 
that are normalizable. Of these six, only three represent 
a different physical form, sn, cn and dn. However, they 
are still related by, 

cn 2 = 1-sn 2 , (Al) 
dn 2 = l-k 2 sn 2 . (A2) 

The six non-normalizable elliptic functions can also be re- 
duced through a phase shift to three with different forms, 
ns, ds, and cs, which are also related by, 

cs 2 = ns 2 - 1 , (A3) 
ds 2 = ns 2 - k 2 . (A4) 

The normalizable and nonnormalizable functions can be 
related through 

sn 2 (iff(l - k 2 ) + z,k) = k 2 ns 2 (z, k) , (A5) 

where K(x) is the complete elliptic integral of the first 
kind. Therefore the square of any elliptic function can 
be related linearly to sn 2 . 

The limits of the sn, cn and dn functions, along with 
the complete elliptic integrals K(k) and E(k) are pre- 
sented in Table [I] The period of the sn 2 , cn 2 and dn 2 
functions is 2K{k). 



APPENDIX B: COMPLETE SOLUTION SET 

It is possible to prove that Eq. Q is a Jacobian elliptic 
integral of the first kind and can therefore be inverted 
to produce the Jacobian elliptic functions. If arbitrary 
parameters are used, Eq. (J7J| becomes, 

f dp = x + x , (Bl) 

J v / A 3 p 3 + A2p 2 +A 1 p + Ao 
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where the Ai's are real constants. The cubic polynomial 
can be factored to give, 



1 



dp = x + x , (B2) 



where at least one of the constants Bi must be real. 
Without loss of generality we may take the real constant 
as B\. The substitution p = y 2 — B\ is then made in 
Eq. (EJ to yield 



Ay = x + x 



^/(y 2 + (B 2 -B 1 ))(y 2 + (B 3 -B 1 )) 



(B3) 

This is the general form of the elliptic integral of the first 
kind 01 an d therefore gives 



Ciel \C 2 y,C 3 ) = x + x 



(B4) 



where the Cj's are constants and el is one of the twelve 
elliptic functions. This can be inverted and p replaced to 
produce 



or finally. 



^/p + B l =C^e\{C^\x + x Q ),C 3 ), (B5) 



p = C 2 2 el 2 (C^(x + x ), C 3 ) - B 1 . (B6) 



Since the square of any elliptic function can be related 
linearly to the square of sn, only one independent solu- 
tion of the form 



p = Aan 2 (bx + x , k) + B , 
need be considered. 



(B7) 
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